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ABSTRACT 
Standard  nonparametric  estimators  of  quantiles  based  on 
order  statistics  can  be  used  not  only  when  the  data  are  i.i.d., 
but  also  when  the  data  are  from  a  stationary,  <j>-mixing  process 
of  continuous  random  variables.   However,  when  the  random  vari- 
ables are  highly  positively  correlated,  sample  sizes  needed  for 
acceptable  precision  in  estimates  of  extreme  quantiles  are  com- 
putationally unmanageable.   A  practical  scheme  is  given,  based 
on  a  maximum  trans formation  in  a  two-way  layout  of  the  data, 
which  reduces  the  sample  size  sufficiently  to  allow  an  experi- 
menter to  obtain  a  point  estimate  of  an  extreme  quantile.   Three 
schemes  are  then  given  which  lead  to  confidence  interval  esti- 
mates for  the  quantile.   One  uses  a  spectral  analysis  of  the 
reduced  sample.   The  other  two,  averaged  group  quantiles  and 
nested  group  quantiles,  are  extensions  of  the  method  of  batched 
means  to  quantile  estimation.   None  of  the  schemes  requires  that 
the  process  being  simulated  is  regenerative. 
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QUANTILE  ESTIMATION  IN  DEPENDENT  SEQUENCES 


1.   INTRODUCTION 

Let   {X  }   be  a  sequence  of  random  variables  which  is 
n 

assumed  to  be  stationary  with  marginal  distribution 

F  (x)  =  P{X   <  x}  .   For  any   p  ,   0  <  p  <  1  ,  the  p-quantile 
j\  n 

of   Fv(x)   satisfies  the  equation   p  =  F  (x  )   if   F  (x)   has 

A  A    p  A 

a  positive  density  in  a  neighborhood  of   x   .   This  paper  is 

concerned  with  the  estimation  of  quantiles   x    from  samples 

X   ,  n=l,...,N.   In  particular  the  paper  is  concerned  with  both 

point  estimates  and  confidence  interval  estimates  for  quantiles 

when  the   X  's   are  highly  correlated.   This  is  the  usual  case 

when  the   X  's   are  the  outputs  from  a  simulation  of  a  stochastic 
n 

system,  for  example  the  waiting  times  in  a  queue. 

In  the  case  of  independent   X  's  ,  methods  for  quantile 

estimation  are  well  known  (e.g.  Connover,  1980,  pp.  71).   Thus 

let   X,  ,   denote  the  order  statistics  from  the  sample  of  size 
(n) 

N  ,  and  let   |_xj   denote  the  greatest  integer  less  than  or  equal 
to   x  .   Then  a  standard  nonparametric  estimator 


x   =  X 


p    A(LpN  +  lj) 


is  a  consistent,  asymptotically  normal,  estimator  of   x    with 
bias  which  is   0(N-1)   and  variance   [ p (1-p) ]/[ f I (x  )N]  +  0(N~2)  , 

A   p 

where   f _,  (x  )   is  the  derivative  of   F„(x)   at   x   .   This  vari- 
X   p  X  p 

ance  can  be  estimated  by  estimating   f  (x  )  ,  but  nonparametric 

confidence  intervals  can  also  be  obtained  for   x    usinq  the 

P 

order  statistics  of  the  sample  (Connover,  1980,  pp.  111-116). 


The  difficulty  with  these  estimation  methods  is  that 
they  require  a  large  amount  of  computer  storage  and  computing 
time  to  sort  the  sample.   Also  for  high  or  low   p  ,  the  order- 
statistic  estimator  might  be  biased  and  very  non-normal.   A 
solution  is  to  use  the  maximum  transformation  (Goodman,  Lewis 
and  Robbins,  1971).   This  not  only  moves  the  problem  back  to  one 
of  estimating  a  median,  but  also  allows  the  use  of  stochastic 
approximation  (Robbins-Monro)  methods.   A  very  satisfactory 
method  of  this  type  using  a  minimal  amount  of  memory  and  giving 
both  point  and  confidence  interval  estimates  for  the  unknown 
quantile   x   has  been  given  by  Robinson  (1975) . 

For  the  case  of  dependent   X  's  ,  the  usual  situation 

^  n 

encountered  in  system  simulation  studies,  quantile  estimation 

is  much  more  difficult  than  in  the  independent  case.   The  point 

estimate   x    is  still  valid;  however  its  variance  is  inflated 
P 

by  a  factor   p(0;x  )  .   Here   p ( 0 ; x  )   is  the  initial  point  on 

the  spectrum  of  the  binary  process   {I  (x  ) }  ,  defined  for  each 

n   to  be   1   if   X   <  x   and  zero  otherwise.   More  directly  it  is 

n  —  p 


p(0;x  )  =  lim  n  Var{(L(x  )  +..+  I  (x  ))/n}  . 
p  1  p         n   p 

n^"°° 


The  nonparametric  confidence  interval  estimation  procedures  are 
no  longer  directly  applicable.   However  it  is  possible  to  esti- 
mate  p(0;x  )   using  methods  of  Heidelberger  and  Welch  (1980, 
1981)  and  to  estimate  the  density   f (x  )   by  standard  methods; 
used  together  these  give  an  estimate  of  the  variance  of   x 

r>  ~o  p 


and  large-sample  confidence  intervals.   The  main  problem  with 
these  order-statistic  estimates  is ,  however ,  that  the  sample 
sizes  required  to  obtain  reasonable  precision  with  the  estimates 
are  prohibitive  when  the   X  's   are  highly  positively  correlated. 

The  basic  scheme  considered  in  this  paper  to  handle  the 
sample  size  problem  for  dependent,  ^-mixing  sequences  is  to  set 
out  the  data  in  a  (conceptual)   v  x  m   array 


Xk,i    Xi  +  (k-l)m  ' 


for   k  =  l,...,v   and   i  =  1,...  m  =  N/v  ,  where   v   is  often 
but  not  necessarily  chosen  so  that   p   =  0.5  .   Then  if   m(and  N) 
is  large  enough  so  that   X  's   which  are   m   apart  are  approxi- 
mately independent,  the  table  can  be  collapsed  by  taking  maxima 
down  the  columns.   The  maxima 


Y .  =   max   X,   .   , 
1    i  ^     k ,  l   ' 
l£k<y 

for   i  =  l,...,m   are  now  a  reduced  sample  whose  q=p   quantile,  y   , 
corresponds  to  the  desired  quantile,   x   ,  of  the   {X  }   process. 
The  main  point  of  this  procedure  is  that  it  gives  a  very  sub- 
stantial sample  size  reduction.   In  addition,  there  is  generally 
slightly  less  overall  correlation  in  the   Y.   sequence  and  this 
counteracts  the  increase  in  variance  which  occurs  with  the  use 
of  the  maximum  transformation  and  order  statistic  estimates  of 

the  quantile   x 

P 


Since  the   Y.'s   (the  maximum  transformed  sample)  are 
still  correlated  several  problems  remain,  notably  estimating 
the  variance  of  the  point  estimate  y   (the  q=pV  sample  quantile  of 

the   Y.'s)   and  obtaining  confidence  interval  estimates  of   x 

1  p 

Three  confidence  interval  schemes  are  considered;  averaged  group 
quantiles,  nested  group  quantiles  and  a  scheme  based  on  esti- 
mating the  probability  density  function  of  the   Y.'s   at   x   , 
and  the  initial  spectral  point,   p(0;y  )   of  the  binary  process 
obtained  from  comparison  of  the   Y.'s   to   y   ,  i.e.   I.  (Ya)  • 

Extensive  empirical  sampling  studies  using   M/G/l   queues 
and  stationary  sequences  of  autocorrelated  exponential  random 
variables  show  that  the  above  schemes  provide  a  reliable  method 
for  estimating  a  quantile  in  a  dependent  sequence   {X  } 


2.   QUANTILE  ESTIMATION  AND  THE  MAXIMUM  TRANSFORMATION 
1 .   Quantile  estimation  for  i.i.d.  sequences. 

Let   X,  ,  .  .  .  ,X  ,...,X    be  a  sample  of  i.i.d.  random 

variables  from  a  continuous  distribution   Fx(x)   with 
probability  density  function   Fx(x)  •   For   0  <  p  <  1   let 


(2.1)  x   =  inf{x  :  Fv(x)  =  p)  =  Fv  (p)   , 

P  A  A 


where   F   (p)   is  the  inverse  of   Fy(x)   with  derivative 

X  A 

l/fv(x  )  .   The  quantity   x  is  the   pth   quantile  of   Fy(x) 
X   p                       P  ** 


r) 


Let   X(1)  <    ,...,<X,,  <  ,...X(N)   be  the  order  stat- 
istics corresponding  to  the  sample.   The  usual  nonparametric 
point  estimator  of   x   is  the   pth   sample  quantile 

ir 


(2.2) 


Xp  -  X  (  |_Np+lJ  )   ' 


0  <  p  <  1 


where   |_zj   denotes  the  integral  part  of   z  .   The  following 
properties  of  x   are  well  known  (David,  1970,  65-67): 


(2.3) 


Bi;  >  =  x  -  f^V   •  EiizPi  +  0 


2fX<V 


N+2 


N2 


(2.4) 


var(x  )  =  a2  =  p(^-p)     +  O 

P  P        Nf2(Xp) 


N 


and 
(2.5) 


N2  (x      -    x    ) 


D 


N(0,    No    ) 


as      N    -* 


Nonparametric  confidence  intervals  for   x    are  based 

P 

on  the  identity  (see  e.g.  Connover,  1980,  pp.  111-116) 


(2.6) 


P{X(L)  i*p  <  X(U)}  =  V 

j—L, 


N 


Pj(i-P)N-j  , 


with  L  <  U   chosen  in  such  a  way  as  to  make  the  probability  as 
close  as  possible  to  the  desired  coverage  (1-a)  . 
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Another  method  for  estimating  quantiles  uses  the  method 
of  stochastic  approximation   (Robbins  and  Monro  (1981) ) .   This 
provides  a  method  of  quantile  estimation  which  avoids  sorting 
and  requires  only  a  minimal  amount  of  storage.   The  asymptotic 
variance  of  the  stochastic  approximation  estimate  is  given  by 
(2.4).   A  very  satisfactory  development  of  this  method  is  given 
in  Robinson  (1975) .   Since  it  is  not  extendable  to  dependent 
data,  it  is  not  described.   The  key  idea  in  the  method,  however, 
is  the  use  of  the  maximum  transformation,  and  since  this  is  of 
use  in  the  present  context,  it  is  described  next. 

2.2.   The  maximum  transformation 

The  maximum  (minimum)  transformation  is  a  computationally 
simple  transformation  and  compaction  of  i.i.d.  data  which  trans- 
forms the  problem  of  estimating  an  extreme  quantile  of  a  random 
variable  into  that  of  estimating  a  more  central  quantile,  e.g. 
the  median.   Thus  assume  that   p  >  1/2  ,  and  let 

Y  =  max (X, , . . . ,X  )  .   (When   p  <  1/2   a  minimum  transformation  is  used) 

Then 

(2.7)        P{Y  <  x  }  =  F  (x  )  =  {F  (x  )}V  =  pV  =  q  . 

—    p        Y    p         A    p 

Thus  the   pth   quantile  of   F    is  the   q  =  pVth   quantile  of  the 

x 

Y  variable.   In  particular  let   v  =  \_in  (1/2)  /in  (p)  \     ;  then   x 
is  approximately  the  median  of  the   Y   variable.   Details  are 
given  in  Goodman,  Lewis,  and  Robbins  (1971). 

There  is  a  price  paid  with  this  transformation  in  infla- 
tion of  the  leading  term  of  the  variance  (2.4).   Thus  let   v 


divide   N  ,  so  that   N/v  =  m  .   Then  the   X   sample  is  reduced 
to  a  sample   Y, , . . . ,Y    by  application  of  the  maximum  transfor- 
mation to  each  successive  non-overlapping  group  of   X.'s   of 

size   v  .   Then  the  density  of  the  Y.'s  is   f  (x  )  =  v  fv(x  ) pV 

i         y   P        •*   P 

and  the  variance  of  the  order  statistic  estimator,   y   , 

of   x    in  the  sample   Y, , . . . ,Y    is 
p  ^      1       m 


v,  ~  ,,       v, 


(2.8,  var,"    )    .    POZE) <Il£L>         .    var(x    )    _ii^ 


NfJ(x  J    v(l-p)pv_1        p   v(l-p)pv  X 

A    p 


The  right-hand  multiplicative  factor  in  this  last  expression  is 
approximately  1.4  if   v   is  chosen  to  make   q  =  p    approximately 
0.5  . 

When  using  stochastic  approximation,  shifting  to  the  neigh- 
borhood of  the   median  is  essential  for  the  method  to  be  a  well- 
behaved  estimation  procedure.   For  estimating  several  quantiles 
simultaneously,  the  method  is  applicable  in  a  nested  scheme 
which  is  very  efficient  with  respect  to  storage  and  speed. 

For   p  <  1/2   a  minimum  transformation  is  used.   Next- 
to-maximum  and  maximum-minimum  schemes  give  more  flexible  and 
robust  schemes  but  are  not  discussed  in  the  present  paper. 

2.3.   Dependent  Data:   General  Considerations  and  Order  Statistics 

Quantile  estimation  in  dependent  data  is  an  order  of 
magnitude  more  difficult  than  in  independent  data.   Fishman 
(1978,  p.   270  )  notes   that  no  satisfactory  solution  exists. 
Iglehart  (1976) ,  Seila  (1976)  and  Moore  (1979)  have  given  special 
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methods  for  processes   {Xn}   with  regenerative  structure,  i.e. 

processes  for  which  there  exist  random  time  points  at  which  the 

process  restarts  probabilistically.   An  example  is  the  waiting 

time  process   {W  }   in  the   M/G/l   queue  which  regenerates  every 

time  a  customer  arrives  to  find  the  queue  empty,  so  that  the 

waiting  time  of  that  customer  is  zero. 

It  is  possible  to  use  the  order  statistic  estimator   x 

P 

given  at  (2.2)  ,  where  we  ignore  for  now  the  problem  of  the 
initial  transient  which  is  usually  encountered  in  simulations. 
Conditions  for  convergence  and  Central  Limit  Theorems  for  sample 
quantiles  are  well  known  (Sen,  1972;  Babu  and  Singh,  1978). 
Thus  let 


(2.9)  In(x)  =1     if     Xn  <  x  , 

=  0     otherwise  ; 


N 
(2.10)  I(x)  =  I       I  (x)/N  , 

n=l 


the  empirical  c.d.f.  at   x  ,  and 


(2.11)  p(0;x)  =  lim  {N  var  I (x) } 


=  I         Yk(x)  , 

k=-°° 

where 

(2.12)  Yk(x)  =  cov{ln(x)'  In+k(x)}         k  =  °'t1"" 

—   m    <    x  <  +  °° 


P{X„    x  ,  Xn+k  <  x}  -  P{Xn  <  x}  P{Xn+k  <  x}  . 


n 


The  notation   p(0;x  )   comes  from  the  fact  that  this 
quantity  is  the  initial  point  (f=0)  of  the  spectrum  of  the 
process   {I  (x  ) }  ; 

1 
(2.13)        P(f;x)  =  I         cos(2Trfk)  yk(x)  ,      -  j   <  f  <  1/2 

k=-°° 


Now  the  Central  Limit  Theorem  for   x_  =  X,,     ,. 

p     (  |.np+lj  ) 

(Sen,  1972)  states  that  if   fv(x)   "*"s  continuous  > 

bounded  and  non-zero  in  some  small  neighborhood  of   x   ,  and  if 

i 

fv(x)   is  bounded  in  this  neighborhood,  then 
x 

N1/2(x   -  x  ) 

(2.14)  £/~ ^ -  N(0,1) 

(P(0;x^))1/Vfy(xn) 
P        A   P 


if  the  process   {X  }   is   ^-mixing  and 


(2.15)  a2   =  p(0;x  )/{f  (x  ) }2 

P         F        F 

is  finite.   For  details  on  ^-mixing  see  Sen  (1972)  and 
Billingsley  (1968);  some  discussion  is  given  in  Section  2.4. 

Regenerative  processes  such  as  the   M/G/l   queuing 
system  waiting  time   {W  }   are  cj) -mixing. 

The  problem  with  using   x    as  an  estimator  in  positively 
correlated  sequences  is  that  the  sample  sizes  required  for 
adequate  precision  are  prohibitively  large.   Both  sorting  times 


10 


and  memory  times  are  then  unrealistic.   A  measure  of  the  infla- 
tion of  sample  size  over  the  independence  case  is  the  ratio  of 
p(0;x  )   to  its  value   p(l-p)   for  the  i.i.d.  process  with 
identical  marginal  distributions.   This  has  been  investigated 
by  Blomqvist  (196  7)  for  the   M/G/l   queue.   For  extreme  quantiles 
of  the   M/M/l   queue  with  traffic  intensity   p  =  0.9   this  ratio 
is  400  .   Greater  ratios  are  possible,  depending  on  the  traffic 
intensity  and  the  skewness  of  the  service  time  distribution. 
Specifically,  the   M/M/l   queue  with   p  =  0.9   requires  a  sample 
size  of  roughly  500,000  customers  to  estimate  the  0.99  quantile 
of  the  waiting  time  distribution  to  within  plus  or  minus  10% 
accuracy.   This  is  derived  from  Table  8  of  Blomqvist  (1967)  ;  more 
precisely  this  is  the  sample  size  required  for  a  90%  confidence 
interval  for  x  „g   to  have  a  relative  half-width  of   0.10  .   For 
the  0.999  quantile  the  required  sample  size  is  approximately 
2,300,000.   Clearly  storing  and  sorting  the  entire  sequence  is 
impractical  in  such  cases.   Actually  to  produce  an  order-statistic 
point  estimate  of   x    requires  storing  only  the  largest   (l-p)N 
values  of  the  sequence.   However,  this  ordering  must  be  dynam- 
ically maintained  as  the  sequence  is  generated,  a  computationally 
expensive  operation.   Additional  storage  is  required  to  estimate 
the  variance  term   p(0;x  ) :   the  positions  in  the  sequence  at 
which  the   (l-p)N   largest  values  occur  must  also  be  saved  in 

order  to  construct  the  binary  sequence   {I  (x  ) }  .   Furthermore, 

J  ^  n   p 

the  Discrete  Fourier  Transform  of  the  extremely  long  sequence 

{I  (x  )  ,  n  =  1,...,N)   must  be  taken  to  accomplish  the  variance 
n   p 

estimation . 

11 


Some  of  the  storage  and  sorting  problems  could  be  re- 
lieved by  using  stochastic  approximation  .   However  properties 
of  the  stochastic  approximation,  particularly  its  asymptotic 
variance,  are  unknown  for  dependent  data.   More  importantly 
direct  application  of  the  maximum  transformation  to  bring  the 
quantile  estimation  down  to  a  problem  of  estimating  the  median 
requires  independence  of  the   xn ' s  • 

We  now  present  a  method  for  using  the  maximum  transfor- 
mation to  achieve  sample  size  reduction  and  a  practical  scheme 
for  point  and  confidence  interval  estimates  with  dependent  data 

2.4.   The  maximum  transformation  in  the  dependent  case 

The  basic  idea  behind  the  use  of  the  maximum  transfor- 
mation is  to  combine  elements  of   X   ,  n  =  1,2,...,  in  a  (f- 
mixing  process  which  are  sufficiently  far  apart  so  as  to  be 
approximately  independent.   To  define  4>-mixing  let   M_    and 

oo 

M      be  respectively  the  o-fields  generated  by   {X.;  i  <  n} 

and    {X.;  1  >  n  +  m}  .   If   E, eM     and   E„cM     ,  then 
1     —  l-oo         2   n+m 

{X  }   is  Q-mixing  if  for  all   n(-°°  <  n  <  °°)   and   m(>  1) 


(2.16)       |P(E2|E1)  -  P(E2)|<   <J>(m)  /       <J>(m)  >  0  , 

where   1  >  <J>(1)  >  (f)(2)  >...  and   lim     cj)  (m)  =0  . 

Thus  we  set  out  the  data  in  a  v^m   array,  where 
m  =  N/v,  as 
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(2.17)      X^.  -  X.+(k_1)m  k=l,...,v;  i=l,...,m  . 

We  assume  for  the  moment  that  elements  in  the   m   columns 

are  independent;  since  we  are  assuming  that   v   is  fixed  and 

the  process   ^xn^   is  (^-mixing  this  will  certainly  be  true  as 

N  (and  m)   get  very  large.   Now  applying  the  maximum  transform 

down  the  columns  we  get  a  reduced  series   Y.  ,  where 

1 


(2.18)        Y.  =   max  X,   .  i  =  1,2,...  ,m 

1    Kk<v  k'1 


Values  of   v   required  to  reduce  the  problem  to  that  of  esti- 
mating the  median  of  the   Y,'s   are  given  in  Table  2.1  . 


Table  2.1 
Size  of   v   for   p   =  1/2 

p    0.900    0.950    0.975    0.980    0.990    0.995    0.998    0.999 
v     6.6     13.5     27.4     34.3     69.0     138.3    346.2    692.8 

We  emphasize  here  that  transformations  to  points  other 

than  the  median  are  possible,  and  are  essential  when  considering 

simultaneous  estimation  of  several  quantiles  or  estimation  of 

the  median  of  the   X  's  .   However  in  what  follows  transforma- 

n 

tion  to  the  median  is  usually  assumed.   There  will  be  an 
attendant  bias  because   p    will  never  quite  be   1/2  . 
The  purposes  of  the  maximum  transformation  are 
(i)   to  reduce  the  attendant  sample  size.   It  should  be  kept  in  mind 
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that  for  higher  quantiles  larger  sample  sizes  are  needed  for 
given  precision.   The  result  is  that  for  fixed  precision  the 
Y.   series  is  roughly  the  same  length  for  all   p  >  0  .   In 
particular  for   p  =  0.999   the  reduction  of  the  sample  size  is 
by  a  factor  of  693  which  is  sufficient  to  reduce  a  series  of 
completely  unmanageable  length  to  one  which  can  clearly  be 
accommodated  in  a  digital  computer  (see  Section  5) . 
(ii)   To  reduce  the  problem  to  estimation  of  a  median  rather 
than  an  extreme  quantile.   Since  stochastic  approximation  is 
not  used  with  the   Y.  series,  this  is  not  essential  for  point 
estimates.   It  is,  however,  helpful  in  obtaining  confidence 
interval  estimates. 

(iii)   To  possibly  reduce  correlation.   It  will  be  shown  that 
the  correlations  in  the   {Y.}   series  are  slightly  less  than 
the  correlations  in  the   {X  }   series.   Though  a  small  effect, 
it  is  usually  sufficient  to  offset  inflation  in  variances  due  to 
use  of  the  maximum  transformation.   To  see  this  assume  that  the 
structure  of   I  (x  )   is  Markovian  with   p(0;x  )  ~  (l+p)/(l-p)  . 
If   p  =  0.99   this  equals  199.   If   p   for  the   Y.   process 
chopped  at  its  median,  i.e.  {I  (xQ  5))  ,  is  reduced  to  only  0.98, 
then  the  new   p(0;x  )  ~  99  . 

sr  p 

Clearly  if   m   is  too  small,  the  assumption  of  inde- 
pendence down  the  column  will  be  violated.   Some  analysis  of 
this  effect  is  possible.   Expanding  the  definition  (2.18)  to 
account  explicitly  for  the  parameter   v  ,   the  number  of   X  's  , 
v   apart,  over  which  the  maximum  is  taken,  so  that 

Y. (v)  =   max   X,   .  ,  i=l,2,...,m  ,  we  have  from  Billingsley 
l<k< v   K  '  x 

(1968,  pp.  174,  20.49)  that 
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(2-19)  |P{Y.(v)  <  xp}  -  pv}|<  vcp(m)  . 

Thus  we  want   m   as  large  as  possible  and,  for  given   v  ,  the 
difference  in  the  probabilities  goes  to  zero  as   m  ->  »  ,  by  the 
definition  of  mixing.   A  slightly  tighter  bound  is  given  by  the 

v-2   , 
Lemma.        |P{Yi(v)  <  x  }  -  pv}  <  <{,  (m)  p  j   pK 

P  k=0 

for   v  >  2  . 


A  straight  forward  proof  of  this  is  obtained  using  inducti 

The  point  estimate  used  for   x    is 

p 


on 


Y~  =  Y 


q     ([mq+lj+l) 


v 


where   q  =  p   .   Since  the   Y_L  '  s   are  dependent  the  properties 
of  this  estimator  are  given  as  at  (2.14)  and  (2.15)  ,  with 

proviso  for  some  peculiarities  in  the  structure  of  the   {Y.} 

i 

process  which  we  discuss  now.   Finite  sample  properties  of  the 
estimator   y    are  studied  by  simulation  in  Section  5. 


2.5.   The  maximum-transformed  process   {Y . }  . 

l 

The  maximum-transformed  process   {Y . }  ,  obtained  from 
the   {X  }   process  by  taking  the  maximum  of   X  's   which  are   m 
apart,  is  not  a  stationary  process  in  the  same  sense  that  the 
{X  }   process  is  stationary.   To  see  this  consider 
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(2.20)     Y   =   max   X     =  max{X1 fXm+1 , . . . ,X     } m+1 }  . 

l<k<v 


It  will  be  correlated  with  successive   Y . ' s  ,  and  this  correla- 
tion can  be  expected  to  decrease  since  the   ^xn^   process  is 
mixing.   However,  the  correlation  will  eventually  increase  because 

(2-21>  Ym  =  max{VX2m Xvm} 

is  coupled  to   Y,   by  the   v-1   dependent  pairs 

{Xm'Xm+l}'---  {X(v-l)m'X(v-l)m+l)   •   Thus  the   {Y.}   process 
is  circular.   However,  it  is  not  stationary  in  a  circular  sense 
because  the  correlation  between   Y,   and 

Y0  =  max{Xn,X  ,„,...,  X,   ,  *  _,_0}   is  through  the   v   dependent 
2         2   m+2       (v-l)m+2 

pairs   {Xl,X2>,  {Xm+1,Xm+2} <x (v-i) ra+l 'X ( v-1) m+2 }  '  andis 

therefore  different  from  the  correlation  between   Y,   and   Y 

1        m 

In  a  circularly  stationary  process  these  correlations  would  be 

the  same.   Towards  the  middle  of  the  process,  say   i  =  m/2  , 

the  lag  one  correlations  will  be  the  same  and  there  will  be  no 

"edge  effect"  because,  if  m   is  large  enough,  the  correlation 

between   Y  /0   and   Y,  ,  and   Y  ,~   and   Y    will  be  zero. 
m/2         1  m/2         m 

Without  belabouring  this  "edge  effect"  and  circularity 
in  the   {Y.}   process  there  are  two  main  points  to  be  made. 

(i)   For  finite   m   it  will  introduce  bias  because,  for 
instance,  methods  using  spectral  techniques  for  estimating 
p(0;x)    in  Section  3  are  based  on  assumptions  of  stationary. 


1G 


Note,  however,  that   {Y.}   is  marginally  stationary,  so  density 
estimation  techniques  for   f  (x)   are  not  affected. 

(ii)   For   m   large  the  effect  will  be  negligible  and 
the  asymptotics  go  through  in  the  usual  way  by  ignoring  a  fixed 
(or  slowly  increasing)  set  of   Y. 's   at  the  beginning  and  end  of 
the  process.   The  usual  assumption  is  that  the  extent  of  this 
set  is  smaller  than  the  extent  of  the  dependence  in  the   {Y. } 
process. 
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3.   CONFIDENCE  INTERVALS  FOR  QUANTILES  IN  DEPENDENT  DATA 

In  the  previous  section  we  gave  a  point  estimator   y 
for  the  quantile  of  the  marginal  distribution  in  a  stationary 
time  series.   In  this  section  we  consider  three  methods  for 
generating  confidence  interval  estimates  for  these  quantiles. 
The  first  method  is  an  extension  of  the  spectral  method  of 
Heidelberger  and  Welch  (1980;  1981).   The  second  and  third  methods 
are  extensions  of  the  method  of  batch  means  (see  e.g.  Law  and 
Carson  (1979)  or  Fishman  (1978)).   These  last  two  methods  also 
give  two  new  point  estimates. 

These  methods  for  generating  confidence  interval  esti- 
mates do  not  rely  upon  use  of  the  max-transform  and  the  case 
v  =  1   corresponds  to  working  with  the  original  time  series. 
Moreover,  if  the  max-transform  is  used  we  do  not  require  that 
p   =  0.50,  i.e.  v  ~  £n.5/£n  p  .   There  are  considerable  compu- 
tational savings,  however,  from  using  some  maximum  transformation. 
This  and  statistical  considerations  will  generally  dictate  a 
maximum  transformation  using  a   v   slightly  less  than  the  median 
transform   v  .   Another  determing  factor  is  that  one  generally 
wants  to  estimate  more  than  one  quantile.   This  factor  will  be 
considered  elsewhere. 

3.1.   A  Spectral  Method 

Let   {Y, ,...,Y  }   be  the  sequence  of  max-transf ormed 

1      m  ^ 

variables  as  defined  by  equation  (2.18).   The  point  estimate  of 

x    is  the   qth   order  statistic  of   {Y,  ,...,Y  }  .  y  ,  where 
P  1      m    Jq 

q  =  p   .   Let   f  (x)   be  the  stationary  density  function  of   Y. 
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and  let  p(0;x)  be  the  stationary  spectral  density  of  the  in- 
dicator function  process   I.(x)   associated  with   {Y.}   at  zero 

i  1 

frequency  as  defined  by  equation  (2.13).   Here  edge  effects 
because  of  the  quasi-circularity  of  the   Y.   process  are  ignored. 

For  large  values  of   m   (see  Sen  (1972)  and  Babu  and  Singh  (1978)) 

^  2       1/2 

^m  (y   -  x  )/(p(0;x  )/f  (x  ) )  /d       has  a  normal  distribution  with 

mean  zero  and  variance  one.   Confidence  intervals  for   x    based 

P 

on  this  Central  Limit  Theorem  are  generated  by  estimating  both 
p(0;x  )  and  f  (x  )  at  the  estimated  quantile  .  This  is  the 
first  method  alluded  to  above. 

The  quantity   p(0;x  )   is  estimated  using  the  spectral 
method  of  Heidelberger  and  Welch  (1980;  1981)  applied  to  the 
sequence   {l.(y  )  ,  i=l,...,m}  .   This  method  uses  least  squares 
to  fit  a  low  order  polynomial  of  degree   d   to  the  logarithm  of 

A 

the  first   K   values  of  the  averaged  periodogram  of   {I. (y  )}  . 

As  suggested  in  Heidelberger  and  Welch  (1980;  1981)  we  used 

K  =  25   and   d  =  2  ,  although  for  extreme  quantiles  in  highly 

congested  queues  with  short  run  lengths,   d  =  3   was  required 

to  produce  valid  confidence  intervals.   This  point  is  discussed 

more  fully  in  Section  4. 

We  found  little  or  no  loss  in  confidence  interval  cov- 

erage  by  using  the  estimated  quantile   y    with  the  sequence 

{I.  (y  )  }   rather  than  using  the  known  quantile   x    with  the 
l  Jq  1  p 

sequence   { I . (x  )}  .   However,  a  proof  of  the  convergence  of 
the  distributional  properties  of  the  periodogram  of   { I ■ ( y  )} 
appears  difficult. 
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The  density   f  (x  )   of  the  max- trans formed  variables 

Y,,...,Y   may  be  estimated  using  a  standard  kernel  density 

1 '      m    2 

estimator  (see  Parzen  (1962)  and  Rosenblatt  (1956)).   Specifically, 
let   W(x)   be  a  weighting  function  such  that 

0  <_   W  (x)  <_  W  <  °°  , 

CO 

(3.1)  /   W(x)dx  =  1  , 


lim   |xW(x) |  =  0 
x  I  -*00 


Let   b(m)   be  a  sequence  of  bandwidth  constants  satisfying 
lim  b  (m)  =  0   and   lim  mb(m)  =  °°  .   For  any   x   the  kernel 

density  estimate  of   f  (x)   is  defined  by 

l   m    l 
(3.2)        f  (x)  =  -   J  T-T-^r   W(  (x-Y.)/b(m)  )  . 
y      m  -=i  b(m)        3 

CO 

1/2 
Under  the  conditions  (3.1)  and  if       4>(n)     <  °°  /  where  the 

n=0 
(J>(n)  '  s   are  defined  at  (2.16)  ,  then   f  (x)   converges  in  prob- 
ability to   f  (x)  . 

In  our  case  we  require  a  density  estimate  at  the  unknown 


/\  ^\ 


point   x   .   Simple  Taylor  series  expansions  show  that   f  (y  ) 
converges  in  probability  to   f  (x  )   if  the  first   k   derivatives, 
W(k) (x)  ,  of   W(x)   satisfy   |w(k) (x) |  <  W  <  «  f  if 

CO 

/  |w(k)(x)|dx       and    lim   |xW(k)(x)|  =  0   and  if 

x  I  ->°° 
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~c 
b(m)  =  Zmm     where   Zm   converges  in  probability  to  a  positive 

constant  and   0  <  c  <  k/2 (k+1)  .   If  the  weighting  function 

W(x)   does  not  satisfy  the  different  ability  conditions  then 

more  delicate  arguments  are  required  to  show  convergence  of 


•\     /\ 


f  (x  )   (see,  for  example,  Robinson  (1975)). 

We  have  experimented  with  two  weighting  functions,  a 
triangular  window 


(l    -    |x|       |x|  <  1  , 


(3.3)  W(x) 

.0  |x|    1  , 


and  a  cosine  window 


(3.4)  W(x) 


x|  <  tt/2  , 

xl  >  tt/2  ; 


Bandwidth  sequences   b (m)   had  the  form   b(m)  =  Z  m 

m 

for  various  powers  of   0  <  c  <  1/2   and  random  variables   Z 

m 

which  ranged  from  constants  to  measures  of  the  spread  of  the 

distribution  of  the   Y.'s  .   We  found  the  density  estimates  to 

3 

be  relatively  insensitive  to  both  the  shape  of  the  weighting 
function  and  the  parameters  of  the  bandwidth  sequence.   For 
small  samples  they  tended  to  slightly  overestimate   f  (x  )   but 
converged  to   f  (x  )   for  large  samples.   In  Section  4  we  report 
the  results  of  experiments  using  the  triangular  window  with 
c  =  1/3   and   Z    equal  to  the  inter-quantile  range  of  the 
Y  .  '  m" 


's  ,  i.e.   b(m)  -  {Y(L.75m+1j)    Y ( L . 25m+l J ) } 
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3.2.   Averaged  and  Nested  Group  Quantiles:   Definition 

We  now  describe  two  methods  of  generating  confidence 
intervals  for  the  quantile   x    which  do  not  require  direct  es- 
timation of   p(0;x  )   for  the   {Y.}   sequence.   These  methods, 
which  we  call  averaged  group  quantiles  (agq)  and  nested  group 
quantiles  (ngq) ,  are  extensions  of  the  method  of  batch  means  to 
quantile  estimation.   Seila  (1976)  considered  a  version  of  agq 
without  the  max  transform   (v=l)   for  regenerative  processes. 

Let   N   be  the  total  sample  size  and  divide   {X, , . . . ,X„} 

1       N 

into   G   non-overlapping  groups  with   mv   observations  in  each 
group   (N  =  Gmv)  .   Define 


(3,5)     X£,j,k    X(£-l)mk+j+(k-l)m     for     £=1,...,G, 

j=l , . . . ,m, 

K —  J.,..., V. 


The  subscript      refers  to  the  group  number,  and  the  data  in 
each  group  can  be  thought  of  as  forming  a  matrix  of  dimension 
v   by   m  .   The  first  row  of  this  matrix  is  formed  from  the 
first   m   observations  in  the  group,  the  second  row  is  made  up 
of  the  second   m   observations  in  the  group,  etc.   For  any  fixed 
values  of  I      and   j   the  points   X   .     and   X   .      (l<k<v) 
are  at  lag   m   (i.e.,  separated  by   m-1   observations)  and  for 

large  values  of   m,  *£  .  x,   X£  j  2"'*,X£  i  k   are  aPProximately 

independent.   Define   Y„ .  =   max  Xn  . 

J         Kk<v  Z>J'k    '    thus   Y£j   1S  the 

jth  max-transformed  variable  in  the   £th   group.   For  large   m 

p{Y£j  £  xp}-3  =  PV   with   lp{Y£j  £  xp}  "  q|  5  <Mm)v   as  in 
Section  2. 

22 


F°r  eaCh   *  =  1 G   let   Y<t(D  i  YM2)---lyMm)   be 

the  order  statistics  of   {Y£1   Yp2  ...  Y   }   and  define 

Yq£  =  Y£(|mq+l|)  '  i,e*   Yq£   is  an  order  statistic  estimate  of 
x    derived  from  the   £th   group.   We  call   y     a  group  quan- 
tile  estimate.   As  in  the  case  of  batch  means,  for  large  values 

of   m  ,   ^yqi  •**'ycfG^   Wl11  be  distributed  as  i.i.d.  normals 

2 
with  mean   x    and  variance   p(0;x  )/mf  (x  )  .   This  suggests 

several  point  and  interval  estimates  for   x 

P 

The  averaged  group  quantile  estimate  is  defined  by 
(3'6)  agq(V  =  k    J=1    Yq£  . 


A  confidence  interval  for   x    is  formed  by  assuming  that 


(3.7)  /G  (agq(xp)  -  xp)/S(agq) 


has  a  students-t  distribution  with   G-l   degrees  of  freedom 
where 

?  1    ^  2 

S  (agq)  =  ^zj      I        (YQi    "  agq(x  ))   . 

i  =  l    M  p 


This  confidence  interval  is  valid  if   {y  ,  . . . ,y  _}   are  i.i.d 

7ql,    ,yqG 

normals  with  mean   x    and  finite  variance. 

P 

The  nested  group  quantile  point  estimate  of   x    is 

defined  to  be  the  median  of   y  ,  .  .  .  ,y  _,  ,  i.e.  if 

Jql  ,     J  qG 

y  ,,>  <  ...  <  y  ,_,  are  the  order  statistics  of   y  ,  . . . ,y    „    , 
Jq(l)  -      -  Jq(G)  ^ql,      J  qG 
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then  the   ngm   is  defined  by  (assuming   G   is  odd) 


(3.8)  ngq(x  )  =  y 


q(L..5G+lJ)  ' 


An  approximate   100*  (1-a) %   confidence  interval  for   x    from 
(2.6)  is 


(3'9)  (Yq(D  ,  yq(U)J  ' 


where  for  any   L  <  U,  a   is  defined  as  at  (2.6) ,  by 

(3.10)  1  -  a  =  I       (£)  (.5)k(.5)G"k  . 

k=L 

This  confidence  interval  is  valid  if   {y  ,  ...,y  _}   are  mutually 

ql ,     JqG  J 

independent  with   P{y  ,  <_   x  }  =  0 . 5  .   For  large  samples  this 
is  guaranteed  by  the  asymptotic  normality  of  order  statistics. 
As   m^-°°,  agq(x  )   and   ngq(x  )   are  asymptotically  un- 
biased and  their  confidence  intervals  are  asymptotically  valid, 
since  the  sequence   {X  }   is  assumed  to  be  mixing.   Notice  that 

if   G  =  1  ,  then  agq(x  )  =  ngq(x  )  =  y    where   y    is  defined 

^  p         P     q  q 

by  equation  (2.20) ,  though  clearly   G  =  1   is  not  reasonable  if 
a  confidence  interval  estimate  is  required.   Furthermore,  as  with 
y      ,  these  point  and  interval  estimates  are  valid  if  the  max- 
transform  is  not  used   (v=l)   in  which  case  one  is  working  with 
group  quantiles  of  the  original  sequence   ^X]}  .   If  the  max- 
transform  is  used  there  is  no  requirement  that   p   =  0.5  . 
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3.3.   Averaged  and  Nested  Group  Quantiles:   Discussi 


on 


There  are  a  number  of  possible  sources  of  error  in  the 

above  two  confidence  interval  procedures.   First,  if   m   is 

small  then,  because  of  dependence,   Hy„  .  <_   x  }   will  differ 

Jo  J  p 

V 

from   q  =  p    by  a  significant  amount  and   E (y  „ )  ^   x   .   Notice 

J  q£  '        p 

that  with   G  >  1  ,   m   for   ngq(x  )   and   agq (x  )   is  smaller 

than  the   m   for   y    and  therefore  this  source  of  error  is 

more  likely  with   ngq(x  )   and   agq(x  )   than  with   y   .   Second, 

p  p  q 

even  if  the   y  ?   sequence  consists  of  i.i.d.  observations,  if 
m   is  small   E (y  . )  ^  x    due  to  the  bias  in  order  statistics  for 
small  samples.   This  bias  will  be  more  severe  at  the  tails  of  the 
distribution.   This  problem  is  again  potentially  more  serious  with 
ngq(x)   and   agq(x  )   than  with   y   .   Third,  the  random  vari- 
ables  y  ,  ...,y     may  be  correlated.   All  three  of  these  fac- 
tors  suggest  that  the  number  of  groups,   G  ,  should  be  as  small 
as  possible.   The  choice   G  =  5   is  the  smallest  value  of   G 

for  which   Ply  ,,,  <  x   <  y  ,_. }  >  0.90   (for   L=1,U=G=5 

7q(l)  -   p    Jq(G) 

the  confidence  level  is  0.9375) .   Given  that   G   should  be  small, 
any  independence  tests  for   y  ,  ...,y     will  have  very  low 
power.   The  sampling  experiments  described  in  Section  4  use 
G  =  5   for  both   agq (x  )   and  ngq (x  )   and  the  correlations 
measured  between  the   y  .'s   have  been  very  low  even  for  strongly 
correlated  sequences  of   X  's   (for  these  experiments  the  sample 
sizes  were  such  that   m  j>   200)  .   Any  residual  correlations 
between  the  group  medians  could  be  further  reduced  by  deleting 
observations  between  groups,  i.e.  by  discarding  observations 
before  starting  the  next  group. 
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3.3.   Theoretical  Comparisons 

It  is  interesting  to  compare  the  expected  confidence 
interval  widths  using  these  three  methods  under  the  assumption 

that  each  method  is  valid.   Let   m  =  N/v   and   m   =  N/vG  =  m/G 

2 

and  assume  (without  loss  of  generality)  that   p(0;x  ) /iruf  (x  )  =  1  . 

Let   G  =  5   and  assume  that   y  ,  . . . ,y  _   are  i.i.d.  normals 

Jql,     ^qG 

2 

with  mean   x    and  variance   p(0;x  )/iruf  (x  )  =  1  .   Then  the 
P  P    <j     P 

expected   ngq   confidence  interval  half-width  is   1.163 
(Pearson  and  Hartley  (1966),  Table  28,  p.  190).   Consider  now 
the  agq  method.   The  t-multiplier  for  a  93.75%  confidence  inter- 
val with  4  degrees  of  freedom  is  2.562  so,  assuming  an  unbiased 
estimate  of  the  standard  deviation,  the  expected  half-width  is 
2.562/V5  =  1.145  .   Finally,  assuming  the  spectral  method  with 

K  =  25   and   d  =  2   produces  an  unbiased  estimate  of  the  point 

2     1/2       — 

estimates'  standard  deviation,   (p(0;x  ) /mf  (x  )  ''     =    1//5  . 

Then  the  confidence  intervals  expected  half-width  is 
2.214//5  =  .99,  since  the  effective  number  of  degrees  of  freedom 
for  the  variance  estimate  was  shown  in  Table  1  of  Heidelberger 
and  Welch  (1981)  to  be  7. 

Thus,  given  the  assumptions,  the  spectral  confidence 
intervals  (with  d  =  2)  will  on  the  average  be  narrower  than  both 
the  agq  and  ngq  intervals  and  the  agq  intervals  will  be  slightly 
narrower  than  the  ngq  intervals.   If  a  cubic  polynomial  is  used 
in  the  spectral  method  then  the  spectral  confidence  intervals  can 
be  shown  to  be  somewhat  wider  than  the  agq  and  ngq  intervals. 
These  relationships  are  empirically  verified  in  Section  4. 
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4.   EMPIRICAL  STUDIES 

All  the  quantile  estimation  methods  given  in  the  previous 
section  are  asymptotically  valid  under  the  4>-mixing  assumption. 
The  crux  of  the  matter  though  is  whether  they  work  with  sample 
sizes  required  for  the  usual  precisions  required  in  systems  simu- 
lations, say  10%  to  5%  as  measured  by  ratio  of  the  confidence 
interval  half-width  to  the  point  estimate.   In  this  section  we 
use  empirical  sampling  (simulation)  to  study  the  bias  and  stand- 
ard deviations  of  the  three  quantile  point  estimators  and  the 
confidence  interval  coverages  and  widths  of  the  three  quantile 
interval  estimators  described  in  Section  3.   The  tests  were  con- 
ducted on  stationary  sequences  of  correlated,  exponentially  dis- 
tributed random  variables  (NEAR(l)  and  GNEAR(l)  processes;  see 
Lawrance  and  Lewis  (1981a)  and  (1981b))  and  on  waiting  time  se- 
quences in  heavily  congested  single  server  queues.   For  each 
process  the  0.50,  0.90,  0.99  and  0.999  quantiles  were  estimated. 

4.1.   Processes  Simulated 

The  NEAR(l)  process   {X  }   with  parameters   a   and   3 
is  defined  as  follows.   Let   {E  }   be  a  sequence  of  i.i.d.  ex- 
ponentials with  mean   1   and  let   {ln}   and   {Knl   be  mutually 
independent  sequences  of  random  variables  for  which 

P{K   =  1}  =  l  -  p{k   =  (l-cOS)  =  <$,  where   6  =  (1-$)/  (1-  (1-a)  3)  / 

n  n 

and      P{ln    =    3>    =    1    -    Pdn    =    0}    =    a    .       Set      XQ    =   EQ    ;    then 


(4.1)  X      =    I    X       ,     +    K    E  ,         n>0 

v         y  n  n    n-1  n   n 
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is  a  Markovian  sequence  of  exponentially  distributed  random 
variables  with  mean  one  and  correlation  structure 


(4.2)     pk  =  corr(Xn,Xn+k)  =  (a3)k  =  pk  ,         k  >  0,  n  >  0 


The  GNEAR(l)  process  yields  a  process  with  exponential 
marginals  and  alternating  positive  and  negative  correlations. 
This  process  is  defined  by   XQ  =  EQ   and 

(4.3)  X   =  K   E   +  I   X'  ,  ,     n>0, 

v/  n     nn     nn-1 

where   X'  =  -  £n(l  -  exp(-  X   ,  ))  .   For  the  GNEAR(l)  process 
n  n-l 

the  lag  one  correlation  is 


(4.4)  p1  =  (aS)r  , 


2 
where   r=l-n/6=-  0.6449   is  the  maximum  negative  cor- 
relation attainable  in  a  bivariate  exponential  distribution. 
Efficient  algorithms  for  generating  the  NEAR ( 1 )  and  GNEAR(l) 
are  given  in  Lawrance  and  Lewis  (1981b). 

The  second  class  of  processes  we  considered  were 
(stationary)  waiting  time  sequences  in   M/G/l   queues,  where 
the  service  time  distribution  has  a  hyperexponential  distribu- 
tion.  Specifically,  let   W   denote  the  waiting  time  of  the 

nth   customer  and  let   {A   ,  n  >  1}   and   {S   ,  n  >  0}   be  the 

n     —  n     — 

i.i.d.  sequences  of  interarrival  times  and  service  times  re- 
spectively.  We  assume  that 
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(4.5) 


P{A   <  x}  =  1  -  exp(-  Ax)  , 


x  >  0  ,   A 


(4.6)     P{Sn  <  x}  =  II^l 


ex 


p(-  y,x) }  +  n2(l  -  exp(-  y2x) }  , 


where   y.   and   II.   are  greater  than  or  eaual  to  zero,  and 
li 


n1   +  n2  =  i 


Let   y  =  (II, /y,  +  n  2/y  2 )     and  let   p  =  A/y  <  1 


be  the  traffic  intensity.   Let   W~   have  the  stationary  waiting 
time  distribution  (see  Kleinrock  (1975) ,  p.  205;  this  distribu- 
tion is  a  probabilistic  mixture  of  an  atom  at   0   and  two  ex- 
ponentials) ;  then  the  waiting  time  sequence   {W  )   is  defined 
by 


(4.7) 


W  _,,  =  (W   +  S   -  A  ,  ,) 

n  +  1      n     n     n  +  1 


n  >  0 


where   x   =  max(x,0) 


We  considered  five  processes  in  testing    the  estimators 


They  are: 

(i) 

(ii) 

(iii) 

(iv) 

(v) 


NEAR(l) 
NEAR(l) 
GNEAR(l) 
M/M/l 


a  =  3  =  0.95; 
a  =  3  =  0.995; 
a  =  3  -  0.995; 
A  =  9,    y,  =  10,  n 


1 ,  p  =  0.90; 


A  =  9,  \i1    =    2,    y2  =  18,  ni  =  0.10,  p  =  0.90; 


M/G/l 

(the  squared  coefficient  of  variation  of  the  service 

time  distribution  is  4.556). 
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The  correlation  structures  of  these  processes  range  from  moderate 
positive  correlation  (NEAR(l),  a  =  3  =  0.95,  pk  ~  . 9  )   to  extreme 
positive  correlation  (NEAR(l),  a  =  3  =  0 . 995 , p  R  ~  . 99   ;  M/M/l 
and   M/G/l)   to  extreme  negative  correlation  (GNEAR(l)).   These 
processes  were  all  simulated  using  the  LLRANDOM  II  random  number 
generating  package  (see  Lewis  and  Uribe  (1981)). 


4.2.   Simulation  Results 

Tables  2-9  report  the  results  of  extensive  simulation 
studies  of  the  quantile  estimation  procedures  detailed  in  earlier 
Sections.   The  notation  used  in  the  tables  is  defined  as  follows. 
For  each  process,  quantile  level   p  ,  value   v   for  the  max- 
transform,  and  run  length  N  ,   R   i.i.d.  replications  were  per- 
formed (R  =  200  for  all  runs,  except  runs  in  which  p  =  0.999, 

and  runs  in  which  p  =  .99  with  large  values  of   m) .   However,  the 
same  random  number  seeds  were  used  for  all  processes,   p,  v  and 

N   so  that  different  entries  in  the  tables  are  not  independent. 

Let   y  (r)  ,  ngq  (x   ,  r)   and   agq(x   ,  r)   denote  the  realiza- 

tions  on  the   rth   replication  of   y   ,  ngq(x  )   and   agq(x  ) 


respectively.   Let   y   ,  ngq   and   agq   denote  the  averages 

Si 

over  the   R   replications  of   Y  (r)  ,  ngq (x  ,r)   and   agq(x  ,r) 
respectively.   For  example 


R  * 
(4.3)  y   =  I      y  (r)/R 

4    r=l 


Let   sd  (y  )  ,  sd(ngq)   and   sd(agq)   denote  the  sample  standard 
si 


deviations  of   y   ,  ngq   and   agq   respectively.   For  example 

Si 
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<4-9)        sd2(y  )  =  I       (y  (r)  -  y  )2/R(R-l)  . 

4     r=l    q        q 

These  empirical  sampling,  point  estimates  and  their  re- 
spective estimated  standard  deviations  may  be  used  to  study  the 
bias  in  the  quantile  estimates  derived  in  this  paper.   For 
example,  from  Table  2,  an  approximate  90%  confidence  interval 
for   E(y  )   in  the  NEAR(l)  process  with   a  =  3  =  0.95,  p  =  0.99, 
N  =  69000   and   v  =  69   is   y   +  1.645  sd(y  )  =  4.599  +  1.645  x  .009 
The  true  value   x  qq  =  4.605  lies  within  this  90%  confidence 
interval,  which  is  (4.584,  4.614). 

(i)   Table  2;  point  estimates  of  quantiles 
Table  2  illustrates  the  effect  of  the  maximum  transfor- 
mation on   y    and  its  standard  deviation,  sd(y  )  .   For  each 

process,  run  length   N   and   p  =  0.90  and  0.99,  y    ancj   sd(v  ) 

q  yq' 

were  computed  using   v  =  1   (no  max  transform)  and   v   chosen 
so  that   pV  ~  0.5  0.   The  run  lengths   N   were  chosen  so  that 
m  =  n/v  =  1000   for  the  NEAR(l)  and  GNEAR ( 1 )  processes  and 
m  =  n/v  =  2000   for  the   M/M/l   and   M/G/l   queues.   These  yield 
conservative  values  of   m.   The  only  case  in  which  the  max- 
transform  introduces  apparent  bias  is  the  M/M/l  queue  with 

p  =  0.90,  v  =  7,  N  =  14000;  in  this  case   y    differs  from   x 

c  J  p  p 

by  2.45  times  the  estimated  standard  deviation.   Since  this  is 
the  maximum  deviation  of  20  experiments,  it  is  probably  not 
significant. 

The  column  labeled  "Actual  sd  inflation"  is  (with   p 

and   N   fixed)  the  ratio   sd(y  |p   ~  0.50)/sd(y  |v=l)  .   This 

j  q  i  c-      ^  j  g  i 

ratio  measures  the  increase  in  standard  deviation  in  dependent 
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sequences  due  to  the  max- trans form.   For  example,  for  the  NEAR(l) 
with   a  =  3  =  0.995  ,  p  =  0.99   and   N  =  69000   this  ratio  is 
.0090/. 0087  =  1.03  . 

The  column  labeled  "Theoretical  sd  inflation  for  i.i.d. 
samples"  is  the  inflation  in  standard  deviation  due  to  the  max- 
transform  for  a  sequence  if  i.i.d.  random  variables;  see 
equation  (2.8).   Recall  that  the  max-transf orm  changes  the  cor- 
relation structure  of  the  process  so  that  its  effect  on  the 
standard  deviation  in  dependent  sequences  may  differ  substantially 
from  its  effect  in  independent  sequences.   In  most  cases  the 
actual  inflation  in  variance  is  very  slight  and  in  several  cases 
a  small  variance  reduction  is  achieved.   Only  for  the   M/M/l 
queue  is  the  inflation  greater  than  what  it  would  be  for  an  i.i.d. 
sequence,  and  even  in  this  case  it  is,  considering  the  storage 
and  computational  savings  of  the  max-transform, an  acceptable 
1.35  . 

From  Table  2  we  conclude  that,  for  sensibly  large  values 
of   m  ,  the  max-transform  introduces  very  little,  if  any,  bias 
and  that  the  inflation  in  variance  is  modest. 

(ii)   Tables  3-8;  confidence  interval  estimates 
Tables  3-8,  in  addition  to  reporting  the  estimated  means 
and  standard  deviations  of  the  quantile  point  estimates, 
compile  the  results  of  confidence  interval  coverage  studies  for 
the  three  quantile  confidence  interval  procedures.   For  each 
process,   p  ,  v  ,  N   and  replication  number,  confidence  inter- 
vals for   x    were  generated  using  the  spectral  method,  nested 
group  quantiles,  and  averaged  group  quantiles  as  described  in 
Section  3. 
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For  these  studies,   G  ,  the  number  of  groups  was  set  at 
5  .   For  this  value  of   G  ,  under  the  assumptions  of  Section  3 
the   ngq   confidence  interval  has  nominal  coverage   0.9375,  i.e. 

P{yq(l)  -  xp    yq(5)}  =  °-9375-   For  comparison  purposes  we  also 
formed  93.75%  confidence  intervals  using  the  spectral  and   agq 
confidence  intervals.   Let   CIr(yg)  ,  CIr(ngq)   and   CIr(agq) 
denote  the  confidence  interval  on  the   rth   replication  using 
the  spectral,  ngq   and   agq   methods  respectively  and  let 


|CIr(y  )|  ,  |CIr(ngq)|   and   |CIr(agq)|   denote  the  widths  of 
these  intervals.   For  each  confidence  interval  procedure 
Tables  3-8  report  the  estimated  average  confidence  interval 
relative  half-width  (labeled  hw/x  )  .   For  example,  for  the 
spectral  confidence  intervals 


R 
(4.10)         hw/x   =  (1/R)  I   I  CI  (y  )|/2x 


These  tables  also  report  the  fraction  of  these  confidence 

intervals  which  actually  contain   x   .   This  fraction  is  called 

P 

an  (estimated)   93.75%   coverage  and  it  should  be  close  to 
0.9375   if  valid  confidence  intervals  are  being  formed.   A  one 
sided  confidence  interval  for  a  coverage,  based  on  the  normal 
approximation  to  the  binomial  distribution,  can  be  used  to  test 
the  hypothesis  that  the  actual  coverage  is  less  than  0.9375  . 
For   R  =  200   (Tables  3,4,5,6  and  part  of  8)  estimated  coverages 
less  than  or  equal  to  0.91  are  significantly  low  at  the  0.90 
level  while  for   R  =  100   (Table  7  and  part  of  8)  estimated 
coverages  less  than  or  equal  to  0.89  are  significantly  low  at 
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the  0.90  level.   However,  in  practice,  the  importance  of  cover- 
age in  judging  the  quality  of  a  confidence  interval  procedure 
depends  very  much  on  the  accuracy  (as  measured  by  relative  half- 
width)  of  the  confidence  intervals.   Thus  low  coverage  is  not  a 
serious  drawback  if  the  procedure  generates  confidence  intervals 
which  are  very  wide  (very  inaccurate)  and  therefore  provide  very 
limited  information  about  the  quantity  being  estimated.   In  fact, 
the  most  useful  information  such  a  confidence  interval  often 
provides  is  that  the  run  length  is  too  short.   As  the  relative 
half-width  decreases,  the  importance  of  coverage  for  the  va- 
lidity of  a  procedure  increases.   We  note  that  it  would  also  be 
possible  to  plot  the  coverage  functions  for  the  spectral  and 
agq  confidence  intervals  (Schruben  (1980));  however,  reporting 
the  individual  93.75%  coverages  is  representative  of  the  cover- 
age function  and  more  compact. 

For  the  spectral  method  and  the   M/M/l   and  M/G/l   queues 
the  coverages  using  both  degrees   d  =  2   and   d  =  3   are  given 
whereas  for  the  NEAR(l)  and  GNEAR(l)  processes  only   d  =  2   is 
given.   For  example,  from  Table  4,  the  coverages  for   x    in  the 
M/M/l   queue  with   p  =  0.90,  v  =  1   and   N  =  14000   for  the 
spectral  method  using  degrees   d  =  2   and   3  ,  nested  group 
quantiles  and  averaged  group  quantiles  are   .905,  .910,  .880 
and   .870   respectively.   The  average  relative  half-widths  are 
.549,  .570,  .429  and  .403   respectively.   Thus  to  increase  the 
precision  from  its  current  value  of  about  50%  to  a  desired  value 
of  10%,  we  estimate  that  the  sample  size  would  have  to  be  in- 
creased by  a  factor  of  25  to  a  total  of  350,000  observations. 
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It  would  clearly  be  impractical  to  store  and  sort  all  of  these 
observations.   Thus  to  increase  the  precision  to  an  acceptable 
level  requires  the  max-transf orm  to  reduce  the  sample  to  a 
manageable  size. 

Notice  that  when  the  coverages  of  all  methods  are  valid 
(for  example  in  Table  5  for  all  run  lengths  of  the  NEAR(l)  and 
GNEAR(l)  processes  and  large  run  lengths  of  the  queueing  processes) 
that  the  relative  widths  of  the  confidence  intervals  are  generally 
as  predicted  in  Section  3,  namely  that 
hw (spectral,  d=2)  <  hw(agq)  <  hw(ngq)  5  hw (spectral,  d=3)  . 

Tables  3  and  4  list  results  of  experiments  for  estimat- 
ing the  median  and  0.90  quantile  respectively  without  using  the 
max-transform  (v  =  1) .   The  run  lengths  used  were  small  enough 
to  accommodate  storing  and  sorting  the  full  sequences. 

In  Table  3  all  coverages  are  either  acceptable  (i.e., 
not  significantly  less  than  0.9375)  or  nearly  acceptable  (i.e., 
wide  confidence  intervals  with  coverage  close  to,  but  signifi- 
cantly less  than  0.9375)  with  the  exception  the  GNEAR(l)  process 
using  the  spectral  method.   The  reason  for  the  low  coverage  in 

A 

this  case  is  that  the  sequence  of  indicator  functions   {I,  (x  ,.)  } 
is  very  nearly  a  deterministic  sequence  of  alternating  zeros 
and  ones  due  to  the  strong  alternating  negative  and  positive 
correlation  structure  of  the  GNEAR(l)  process.   This  phenomenon 
appears  to  be  peculiar  to  the  median  (the  coverage  for  the  0.90 
quantile  using  the  spectral  method  in  Table  4  is  acceptable)  but 
should  be  kept  in  mind  when  dealing  with  processes  having  this 
type  of  correlation  structure. 
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The  coverages  for  the  0.90  quantiles  in  Table  4  are  all 


acceptable  or  nearly  acceptable.   However,  agq   and   ngq   exhibit 
substantial  bias  in  the  M/G/l   queue.   In  both  Tables  3  and  4 
notice  the  very  large  mean  relative  half-widths  in  the   M/M/l 
and  M/G/l   queues,  again  implying  the  need  for  the  max-transform 
to  deal  with  the  much  larger  sample  sizes  required  for  estimates 
of  acceptable  precision. 

Tables  5,  6  and  7  compile  results  of  experiments  for 
estimating  the  0.90,  0.99  and  0.999  quantiles,  respectively, 
using  the  max-transform  with   v   chosen  so  that   p   ~  0.50  . 

The   ngq   and   agq   confidence  interval  estimate  coverages 
are  generally  acceptable  throughout  these  tables.   However,  as 


in  the  case  of   v  =  1  ,   ngq   and   agq   show  that  there  is  sub- 
stantial bias  in  the   ngq   and   agq   point  estimates  in  the 
M/G/l   queue  for  small   N.  For   N  =  224,000   this  bias  has  dis- 
appeared even  though  the  precision  (about  20%)  is  still  relatively 
high.   The  estimates   x    show  little,  if  any,  bias.   The  cover- 
ages of  the  confidence  interval  estimates  in  the  M/G/l   queue 
are  high  despite  the  bias  due  to  the  extreme  width  of  the  confi- 
dence intervals  at  the  small  run  lengths. 

For  the  spectral  method  the  coverages  for  the  NEAR(l) 
and  GNEAR(l)  processes  are  acceptable.   However,  the  small  sample 
coverages  for  the  queues,  particularly  the   M/G/l   queue,  are 
disappointing.   Even  using  a  cubic  polynomial  to  fit  the  log  of 
the  periodogram  does  not  yield  acceptable  coverages  for  the 
M/G/l   queue  when   N   is  small.   The  low  coverage  is  explained 
by  the  fact  that  the  spectrum  of  the  max-transformed  variables 
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{Ykl   is  still  very  steep  near   f  =  0   (i.e.,  the   Y.'s   are 
highly  correlated)  even  for  large  values  of   v  .   Since   N  =  mv  , 
if   v   is  large  then   m   is  small  and  the  least  squares  fits  to 
the  log  of  the  periodogram  is  done  over  a  relatively  large  fre- 
quency range.   If  log  (p(f  ;  x  )  )  is  very  steep  then  it  may  not 
be  well  approximated  by  a  low  order  polynomial  over  a  wide  fre- 
quency range. 

Figure  1  illustrates  this  problem  in  the  case  of  the 
NEAR(l)  process  with   a  =  3  =  0.995  .   The  log  of  the  spectrum 
of  the  max-transformed  variables  (v  =  693)  is  seen  to  be  still 
quite  steep  near   f  =  0  .   Compare  this  to  the  almost  flat  log 
of   PD(f)  ,    the  spectrum  of  the  process  of  batch  means   X^(k)  , 

B  B 

where 

,      kB 
(4.11)  X  (k)  =  ±      I  X 

B       B  j=(k-l)B+l   : 

These  batch  means  would  be  used  to  place  a  confidence  interval 
on   E  (X  .  )  (see  Heidelberger  and  Welch  (1981)). 

The  problem  can  be  avoided  by  making   v   smaller,  and 
thus   m   larger.   The  least  squares  fit  is  then  done  over  a 
narrower  frequency  range  and  the  coverage  improves.   This  effect 
is  seen  by  comparing  the  coverages  for   x  Qfi   in  the   M/G/l 
queue  with   N  =  14000   and   v  =  1  (Table  4)  to  those  using  v  =  7 
and   N  =  14000   (Table  5) .   (Again  note  that  the  precision  is 
unacceptably  high) .   It  is  also  seen  in  Table  8  which,  for  the 
M/M/l   and   M/G/l   queues,  compares  the  coverages  for   x  gg 
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Effect  of  the  maximum  transform  and  batching 
on  the  log-spectrum.   Figure  la  gives  the  log 
of  the  spectrum  of  a  NEAR(l)  process  with 
a  =  3  =  0.995.   Figure  lb  shews  the  log  of 
the  spectrum  after  the  maximum  transformation 
with  v  =  693.   A  small  but  significant  effect 
can  be  seen 'at  f  =  0.   Figure  lc  shows  the 
log  of  the  spectrum  after  batching,  with  batch 
sizes  of  v  =  693.   Batching  reduces  the  initial 
point  of  the  spectrum  much  more  than  the  maximum 
transformation. 
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Figure  1.   Effect  of  the  maximum  transform  and  batching 
on  the  log-spectrum.   Figure  la  gives  the  log 
of  the  spectrum  of  a  NEAR(l)  process  with 
a   =    p  =  0.995.   Figure  lb  shews  the  log  of 
the  spectrum  after  the  maximum  transformation 
with  v  =  693.   A  small  but  significant  effect 
can  be  seen  at  f  =  0.   Figure  lc  shows  the 
log  of  the  spectrum  after  batching,  with  batch 
sizes  of  v  =  693.   Batching  reduces  the  initial 
point  of  the  spectrum  much  more  than  the  maximum 
transformation. 
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using   v  =  69   to  those  using   v  =  17   and   v  =  8  ,  and  the 
coverages  for   x  g9g   using   v  =  693   to  those  using   v  =  173 
and   v  =  86   (the  slight  differences  in   N   are  to  make   m 
highly  composite  which  increases  the  computational  efficiency  of 
the  Discrete  Fourier  Transform  used  to  calculate  the  periodogram) . 
Thus  the  full  storage  and  computational  savings  of  the  max- 
transform  (i.e.  moving  the  quantile  back  to  the  median)  cannot 
be  realized  in  highly  correlated  sequences  without  sacrificing 
confidence  interval  coverage.   Significant  savings  can  be  realized 
however.   For  example  the  max-transf orm  with   v  =  173   reduces 
the  1,384,000  observations  to  a  sample  size  of  just  8000  for 
estimating   x  _„«   in  the   M/G/l   queue. 

In  Table  8  two  simulations  are  also  shown  with 
N  =  5,536,000   and   v  =  173   for  both  the   M/M/l   and   M/G/l 
queues  for  the  extreme  case  of  the  .999  quantiles.   This  sample 
size  is  large  enough  to  obtain  almost  10%  precision.   To  within 
the  precision  of  the  empirical  sampling  experiment  (R  =  100)  ,  it 
should  be  noted  that  all  point  and  confidence  interval  estimates 
are  valid.   Note  that  the  maximum  transformation  reduces  the 
sample  size  from  5,536,000  to  32,000! I 
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5.   SUMMARY  AND  CONCLUSIONS 

In  this  paper  we  have  investigated  point  and  interval 
estimates  for  a  single  quantile  in  a  fixed  length  sequence  of 
dependent  observations.   For  large  sample  sizes  it  becomes  im- 
practical to  store  and  sort  the  entire  sequence.   For  extreme 
quantiles,  these  limitations  can  be  overcome  by  using  the  maxi- 
mum transformation  which  requires  storing  only  a  sequence  of 
maxima.   This  sequence  is  defined  by  laying  the  data  out  into  a 
v   by   m   array  and  storing  only  the  maximum  element  in  each 
column.   Storage  requirements  are  thus  reduced  by  a  factor  of   v  . 
Observations  at  lag   m   are  assumed  to  be  independent  and  the 
max-transform  changes  the  problem  of  estimating  the   p   quantile 
of  the  original  sequence  into  one  of  estimating  the   p    quantile 
of  the  max-transf ormed  sequence. 

Three  confidence  interval  methods  which  can  exploit  the 
sample  size  reduction  produced  by  the  max-transf ormation  were 
described  and  tested;  the  spectral  method  and  two  extensions  to 
the  method  of  batch  means  called  nested  group  quantiles  (ngq) 
and  averaged  group  quantiles  (agq) .   The  three  methods  produced 
valid  confidence  intervals  if  the  sample  sizes  were  large  enough 
to  produce  10%  precision  in  the  estimates. 

Problems  which  were  not  addressed  in  this  paper  and  are 
the  subject  of  ongoing  research  involve  questions  of  multiple 
quantile  estimation  and  the  incorporation  of  these  fixed  run 
length  procedures  into  sequential  procedures  for,  say,  simulation 
run  length  control.   Among  the  issues  involved  here  are: 
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(i)   Testing  the  adequacy  of  the  spacing  parameter   m 
so  that   X,   and   X,     are  approximately  independent. 

(ii)   Organizing  the  data  so  that   m   can  be  increased 
as  the  sample  size  increases. 

(iii)   Organizing  the  data  so  that  multiple  quantiles  can 
be  efficiently  estimated,  for  example,  estimating  a  0.50  and 
0.999  quantile  for  the  same  set  of  observations  (this  probably 
involves  the  use  of  a  max-min  transformation) . 

(iv)   Determining  the  sensitivity  of  these  schemes  to 
an  initial  transient  which  is  often  encountered  in  simulation. 
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Table  2.  Comparison  of  point  estimates  and  standard  deviations  of  point  estimates 

using  full  sample  order  statistic  and  max- trans formed  sample  order  statistic 
(R  =  200  replications) .   Here   q  =  pv  .   The  fact  that  the  actual  sd  inflation 
is  generally  less  than  the  theoretical  sd  inflation  for  i.i.d.  sequences  is 
because  of  the  decrease  in  correlation  in  the  series  induced  by  the  maximum 
transformation. 

Process 

P 

x 

P 

N 

V 

\ 

sd(yg) 

Actual 
sd 
inflation 

Theoretical  sd 
inflation  for 
i.i.  d.  sequences 

NEAR(l) 
a=3=0.95 

0.90 
0.90 

0.99 
0.99 

2.303 
2.303 

4.605 
4.605 

7,000 
7,000 

69P00 
69,000 

1 

7 

1 
69 

2.304 
2.320 

4.610 
4.599 

.0104 
.0106 

.0087 
.0090 

1.02 
1.03 

1.18 
1.20 

NEAR(l) 
a= 3=0. 995 

0.90 
0.90 

0.99 
0.99 

2.303 
2.303 

4.605 
4.605 

7,000 
7,000 

69,000 
69,000 

1 

7 

1 
69 

2.344 
2.340 

4.540 
4.607 

.0325 
.0356 

.0275 
.0299 

1.10 
1.09 

1.18 
1.20 

GNEAR(l) 

a=3=0.995 
negative 
correlation 

0.90 
0.90 

0.99 
0.99 

2.303 
2.303 

4.605 
4.605 

7,000 
7,000 

69^00 
69,000 

1 

7 

1 
69 

2.282 
2.291 

4.605 
4.600 

.0132 
.0121 

.0046 
.0050 

0.92 
1.09 

1.18 
1.20 

M/M/1 
p=0.90 

0.90 
0.90 

0.99 
0.99 

2.197 
2.197 

4.500 
4.500 

14,000 
14P00 

138p00 
133000 

1 

7 

1 
69 

2.222 
2.289 

4.514 
4.558 

.0277 
.0375 

.0368 

.0455 

1.35 
1.24 

1.18 
1.20 

M/G/l 
P=0.90 

0.90 
0.90 

0.99 
0.99 

6.311 
6.311 

13.130 
13.130 

14,000 
14P00 

138P00 
138,000 

1 

7 

1 
69 

6.466 
6.572 

13.182 
13.398 

.1682 
.1759 

.2077 
.1975 

1.05 
0.95 

1.18 
1.20 
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Table  3. 

i 

Point  estimates 
coverages,  and 

three  methods, 
of  replications 

;,  standard  c 
average  conJ 
Here  p  =  0. 
;,  is  200. 

leviations  of  point  estimates, 
iidence  interval  half-width/x 
50,  v  =  1,  G  =  5,  and  R,  theF 

93.75% 
for  the 
number 

e 
icess 

X 

P 

N 

Spectral  Method 

Nested  Group 
Quantiles 

Averaged  Group 
Quantiles 

X 

P 

sd(x  ) 
P 

degree  of  polynomial 
3          2 

ngq 
sd(ngq) 

coverage 

hw/x 
P 

ngq 
sd(agq) 

coverage 

hw/x 
P 

coverage 

hw/x 
P 

coverage 

hw/x 
P 

5AR(1) 

3=0.95 

0.693 

7p00 

.696 
(.003) 

.935 
(.149) 

.700 
(.004) 

.935 
(.180) 

.706 
(.003) 

.945 
(.167) 

i 

EAR(l) 
3=0.995 

0.693 

7,000 

.735 
(.012) 

.865 
(.517) 

.744 
(.015) 

.905 
(.598) 

.809 
(.013) 

.910 
(.565) 

EAR(l) 
13=0.995 
jative 
irelation) 



0.693 

7,000 

.692 
(.001) 

.795 
(.006) 

.693 
(.002) 

.945 
(.112) 

.689 
(.002) 

.985 
(.102) 

/M/l 
0.90 

0.588 

14P00 

.601 
(.006) 

.950 
(.452) 

.950 
(.352) 

.598 
(.007) 

.945 
(.443) 

.640 
(.008) 

.940 
(.413) 

I/G/1 
=  0.90 

1.545 

14P00 

1.607 
(.033) 

.950 
(.878) 

.940 
(.647) 

1.562 
(.036) 

.930 
(1.013) 

1.943 
(.055) 

.940 
(.961) 

1 
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Table  4.   Point  estimates,  standard  deviations  of  point  estimates,  93.75% 

>r  thi 

v  -  1, 


and  average  confidence  interval  half-width/x  for  the 


three  methods,  no  maximum  transform.   Here  p  -  0.90, 
G  =  5  and  R,  the  number  of  replications,  is  200. 


Process 

i 

X 

P 

N 

Spectral  Meth 

3d 

Nested  Group 
Quantiles 

Averaged  Group 
Quantiles 

X 

P 

sd(x  ) 
P 

degree  of  p< 
3 

Dlynomial 
2 

ngq 
sd(ngq) 

coverage 

hw/x 
P 

ngq 
sd(agq) 

coverage 

hw/x 
P 

coverage 

hw/x 
P 

cove rag 

hw/x 
P 

NEAR(l) 
a=3=0.95 

2.303 

7,000 

2.304 
(.010) 

.930 
(.129) 

2.287 
(.013) 

.925 
(.156) 

2.296 
(.010) 

.915 
(.145) 

NEAR(l) 
a= 3=0. 995 

2.303 

7p00 

2.344 
(.032) 

.935 
(.468) 

2.189 
(.038) 

.955 
(.469) 

2.264 
(.030) 

.940 
(.435) 

GNEAR(l) 
a= 3=0. 995 
(negative 
correlation) 

2.303 

7,000 

2.305 
(.010) 

.925 
(.141) 

2.282 
(.013) 

.955 
(.181) 

2.267 
(.011) 

.955 
(.168) 

M/M/1 
p=0.90 

2.197 

14,000 

2.222 
(.028) 

.910 
(.570) 

.905 
(.549) 

2.066 
(.027) 

.880 
(.429) 

2.229 
(.027) 

.870 
(.403) 

M/G/1 
p=0.90 

6.311 

14P00 

6.466 
(.168) 

.890 
(1.069) 

.845 
(.669) 

5.266 
(.090) 

.865 
(.496) 

5.843 
(.099) 

.815 
(.466) 
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Table  5. 

Point  estimates 

,  standard 

deviations  c 

3f  point 

estimates , 

93.75% 

coverages,  and  < 

average  con 

fidence  interval  half-width/x 

for  the 

three  methods. 

As  in  Table  4,  p  ■■     .90,  but  here  v  =  7; 

thus  pv  ' 

:-  0.50  . 

Also  G 

=  5  and  ] 

\,    the  numb 

er  of  replications,  is  200. 

Spectral  Method 

Nested  Group 

Average 

d  Group 

Quantiles 

Quantiles 

1 

^q 

degree  of 
3 

polynomial 
2 

ngq 

coverage 

ngq 

| 

coverage 

coverage 

coverage  : 

Ip  cess 

X 

P 

N 

sd(y  ) 

q 

hw/x 
P 

hw/x 
P 

sd(ngq) 

hw/x 
P 

sd(agq) 

hw/x 

P 

i 

|ear(1) 

2.303 

7,000 

2.320 

.945 

2.310 

.940 

2.328 

.940 

13=0.95 

i.011) 

(.156) 

(.025) 

(.270) 

(.012) 

(.164) 

28,000 

2.304 

.940 

2.30  3 

.920 

2.305 

! 

.915 

(.006) 

(.075) 

(.007) 

(.089) 

(.006) 

(.083) 

112,000 

2.300 

.930 

2.296 

.920 

2.300 

.9  35 

(.003) 

(.037) 

(.005) 

(.044) 

(.003) 

(.040) 

EAR(l) 

2.303 

7,000 

2.340 

.800 

2.253 

.940 

2.414 

.945   ; 

3=0.995 

(.036) 

( . 342) 

(.038) 

(.519) 

(.033) 

(.488) 

28,000 

2.311 

.935 

2.297 

.960 

2.350 

.940 

(.017) 

(.246) 

(.020) 

(.270) 

(.027) 

(.200) 

112,000 

2.289 

.970 

2.301 

.945 

2.309 

.945 

(.008) 

(.127) 

(.010) 

(.242) 

(.008) 

(.752) 

KEAR(l) 

2.303 

7,000 

2.291 

.960 

2.316 

.935     2.287 

.940 

=3=0.995 

(.012) 

(.153) 

(.013) 

(.192) 

(.012) 

(.178) 

28,000 

2.306 

.940 

2.312 

.940 

2.299 

.935 

egative 

(.006) 

(.078) 

(.007) 

(.007) 

(.000) 

(.000) 

orrelation) 

112,000 

2.304 

.925 

2.306 

.925 

2.304 

.935 

(.003) 

(.040) 

(.003) 

(.047) 

(.003) 

(.043) 

M/M/l 

2.197 

14,000 

2.289 

.940 

.815 

2.206 

.955 

2.331 

.930 

=0.90 

(.037) 

(.500) 

(.306) 

(.055) 

(.452) 

(.052) 

(.422) 

56,000 

2.231 

.975 

.950 

2.220 

.945 

2.316 

.935 

(.016) 

(.325) 

(.240) 

(.020) 

(.296) 

(.020) 

(.274) 

224,000 

2.206 

.965 

.950 

2.204 

.955 

2.227 

.945 

(.008) 

(.142) 

(.176) 

(.009) 

(.255) 

(.008) 

(.127) 

M/G/l 

6.311 

14,000 

b.512 

.695 

.525 

5.657 

.865 

6.112 

.860 

=0.90 

(.176) 

(.399) 

(.219) 

(.725) 

(.500) 

(.220) 

(.472) 

56,000 

6.373 

.950 

.870 

6.184 

.940 

6.680 

.935 

(.079) 

(.471) 

(.500) 

(.079) 

(.42-;) 

(.085) 

(.417) 

224,000 

6.388 

.970 

.955 

6.367 

.965 

6.559 

.960 

(.042) 

(.260) 

(.209) 

(.047) 

(.245) 

(.047) 

(.227) 

Table  6 

Point  estimates,  standard 
coverages,  and  average  con 
three  methods.   Here  p  =  0 
q  =  pv  ~  0.50  .   Also  G  = 
is  200. ~  However  R  =  100 

deviations  of  point  estimates,  93.75% 
fidence  interval  half-width/x  for  the 
.99,  v  =  69,  so  that         P 
5  and  the  number  of  replications,   R, 
if  the  run  is  marked   t  . 

Spectral  Method 

Nested  Group 
Quantiles 

Averaged  Grouj 
Quantiles 

"yq 

degree  of 
3 

polynomial 
2 

ngq 

coverage 

ngq 

coven 

coverage 

coverage 

Process 

X 

P 

N 

sd(9q) 

hw/x 
P 

hw/x 
P 

sd(ngq) 

hw/x 
P 

sd(agq) 

hw/x 

I 

NEAR(l) 

a=3=0.95 

4.605 

69,000 

4.599 

{.009) 

.960 
{.063) 

4.606 

{.011) 

.920 
{.073) 

4.627 

(.000) 

.91! 

{.06; 

138,000 

4.608 
{.006) 

.945 
{.043) 

4.604 
{.007) 

.950 
(.054) 

4.608 
(.000) 

.95C 

(.056 

276,000 

4.604 
{.004) 

.930 
{.031) 

4.603 
(.005) 

.945 
(.037) 

4.608 

(.005) 

.95C 

(.034 

NEAR(l) 

4.605 

69,000 

4.607 
{.030) 

.890 

{.180) 

4.606 
(.052) 

.965 
{.248) 

4.692 
{.029) 

.965 
(.22£ 

138,000 

4.653 
{.023) 

.930 
{.140) 

4.653 

{.027) 

.905 
(.273) 

4.685 
{.024) 

.925 
(.206 

2  76,000 

4.620 

{•014) 

.940 
{.099) 

4.613 
(.020) 

.945 

{.111) 

4.629 

(.024) 

.935 

{.10k 

GNEAR(l) 
a= 3=0. 995 

4.605 

69,000 

4.600 

{.005) 

.955 
{.035) 

4.596 
(.000) 

.950 
{.040) 

4.594 
(.005) 

.94C 

(.03? 

(negative 

138,000 

4.605 
{.004) 

.920 
{.023) 

4.605 
(.005) 

.925 
(.025) 

4.601 
(.004) 

.91C 
{.026 

correlation) 
1 

276,000 

4.608 
{.003) 

.950 

{.017) 

4.604 
(.005) 

.930 
(.020) 

4.605 

(.003) 

.935 
(.02fl 

M/M/l 

4.500 

138,000 

4.558 
{.045) 

.865 
{.262) 

.730 
(.253) 

4.351 
{.035) 

.920 
(.207) 

4.538 
(.035) 

.905 
{.250 

276,000 

4.513 

{.027) 

.955 
{.236) 

.895 
{.148) 

{.032) 

.930 
{.222) 

4.598 
(.025) 

.930 

{.207 

552,000 

4.525 
{.020) 

.955 
{.190) 

.945 
(.132) 

4.511 
{.024) 

.945 
(.275) 

4.608 
{.023) 

.935 

{.166 

t 

2,208,000 

4.481 
{.012) 

.980 

{.084) 

.960 
(.057) 

4.495 
(.025) 

.960 
(.070) 

4.510 

(.023) 

.950 
(.075 

M/G/l 

13.130 

138,000 

13.398 
{.198) 

.675 
{.243) 

.490 
{.129) 

11.713 
(.225) 

.820 
{.288) 

12.174 
{.119) 

.815 
(.270 

2  76,000 

13.359 
{.154) 

.845 
{.270) 

.635 
(.155) 

12.377 
{.107) 

.920 
(.257) 

12.934 
(.200) 

.900 
(.207 

552,000 

13.271 
{.106) 

.935 
{.276) 

.815 
{.164) 

12.927 
(.203) 

.910 
{.265) 

13.480 

(.20i) 

.915 
(.247 

t 

2,208,000 

13.102 
{.062) 

.960 

{.163) 

.930 
(.i23) 

13.010 

{.075) 

.920 
(.233) 

13.241 
(.000) 

.900 
{.125 

Table  7.   Point  estimates,  standard  deviations  of  point  estimates,  93.75% 
coverages,  and  average  confidence  interval  half-width/x   for  the 
three  methods.   Here  p  =  0.999,  v  =  693,  so  that       P 
q  =  PV  ^  0.50  .   Also  G  =  5  and  the  number  of  replications,   R, 
is  100. 


1 

X 

P 

N 

Sp 

lectral   Method 

Nested    Group 
Quantiles 

Averaged   Group 
Quantiles 

'  ':! 

sd(yq) 

degree   of   pc 
3 

)lynomial 
2 

ngq 

sd(ngq) 

coverage 

hw/x 
P 

ngq 
sd(agq) 

coverage 

hw/x 
P 

Ipcess 

coverage 

hw/x 
P 

coverage 

hw/x 
P 

JEAR(l) 
13=0.95 

6.908 

693P00 

6.923 
(.011) 

.940 
(.034) 

6.908 
(.012) 

.940 
(.038) 

6.912 
(.010) 

.940 
(.035) 

EAR(l) 
3=0.995 

6.908 

693P00 

6.950 
(.033) 

.960 
(.117) 

6.937 
(.040) 

.920 
(.138) 

7.018 
(.037) 

.920 
(.128) 

|ear(d 

=  3=0.995 

»ative 

^relation) 

6.908 

693£00 

6.906 
(.005) 

.930 
(.015) 

6.899 
(.006) 

.930 
(.018) 

6.899 
(.005) 

.930 
(.016) 

, 

1/M/l 
=0.90 

6.800 

1386P00 

6.835 
(.058) 

.900 
(.181) 

.730 
(.108) 

6.648 
(.046) 

.930 
(.172) 

6.790 
(.046) 

.940 
(.160) 

1/G/l 

=  0.90 

19.94 

1^  86,000 

20.693 
(.324) 

.700 
(.161) 

.520 
(.093) 

18.589 
(.171) 

.890 
(.211) 

19.063 
(.157) 

.840 
(.198) 
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Table  8 

M/M/1  queue.   Point  estimates,  standard  deviations  of  point 
estimates,  93.75%  coverages,  and  average  confidence  interval 
half-widths/x  comparing  q  =  pv  ~  0.50,  q  =  pv  ~  0.84  and 
q  =  pv  ~  0.92   for  different  sample  sizes.   Here  p  =  0.99  and 
p  =  0.999  and  G  =  5.   The  number  of  replications  is   R  =  200 
for  runs  marked  *    ,  and  R  =  100  for  runs  marked  t  . 

P 

Spectral  Method 

Nested  Group 
Quantiles 

Averaged  Group 
Quantiles 

^q 

degree  of 
3 

Dolynomial 
2 

ngq 

coverage 

ngq 

coverag 

coverage 

coverage 

Process 

X 

P 

N 

sd(y  ) 

hw/x 
P 

hw/x 
P 

sd(ngq) 

hw/x 
P 

sd(agq) 

hw/x 

P 

M/M/1    * 
p=0.90 

0.99 
A. 500 

138,000 
69 

4.558 
(.045) 

.865 
(.262) 

.730 
(.153) 

4.351 
(.035) 

.920 
(.267) 

4.538 
(.035) 

.905 
(.250) 

* 

0.99 
4.500 

136,000 
17 

4.546 

(.040) 

.925 
(.375) 

.880 

(.246) 

4.258 
(.029) 

.895 
(.245) 

4.409 
(.029) 

.890 
(.222) 

* 

0.99 
4.500 

128,000 
8 

4.596 
(.042) 

.900 
(.384) 

.870 
(.306) 

4.317 
(.035) 

.895 
(.260) 

4.470 
(.557) 

.875 
(.245) 

t 

0.999 
6.800 

1,386,000 
693 

6.835 

(.058) 

.900 

(.181) 

.730 

(.108) 

6.648 

(.046) 

.930 
(.272) 

6.790 
(.046) 

.940 

(.160) 

t 

0.999 
6.800 

1,384,000 
173 

6.855 
(.059) 

.910 
(.284) 

.870 
(.178) 

6.578 
(.055) 

.930 

(.275) 

6.732 
(.555) 

.910 
(.255) 

t 

0.999 
6.800 

1,376,000 
86 

6.856 
(.056) 

.940 
(.232) 

.890 

(.181) 

6.657 
(.054) 

.900 
(.156) 

6.754 
(.545) 

.880 

(.144) 

t 

0.999 
6.800 

5,536,000 
173 

6.797 
(.024) 

.970 
(.112) 

.950 

(.088) 

6.743 
(.027) 

.930 
(.256) 

6.833 
(.028) 

.940 

(.255) 

M/G/1     * 
p=0.90 

0.99 
13.13C 

138,000 
69 

13.398 
(.198) 

.675 
(.243) 

.490 
(.129) 

11.713 
(.225) 

.820 
(.288) 

12.174 

(.225) 

.815 

(.275) 

* 

0.99 
13.13C 

136,000 
17 

13.240 

(.202) 

.850 
(.462) 

.725 
(.274) 

11.467 
(.125) 

.795 
(.254) 

12.011 
(.124) 

.755 
(.277) 

* 

0.99 
13.13C 

128,000 
8 

13.011 

(.200) 

.805 
(.758) 

.730 
(.489) 

11.336 
(.120) 

.755 
(.265) 

11.762 
(.135) 

.720 
(.266) 

t 

0.999 
19.94 

1,386,000 
693 

20.693 
(.324) 

.700 
(.161) 

.520 
(.093) 

18.589 
(.171) 

.890 
(.211) 

19.063 

(.257) 

.840 

(.198) 

t 

0.999 
19.94 

1,384,000 
173 

19.951 
(.230) 

.890 
(.294) 

.800 

(.170) 

18.188 
(.159) 

.880 
(.198) 

18.638 

(.137) 

.880 

(.185) 

t 

0.999 
19.94 

1,376,000 
86 

19.844 
(.239) 

.820 
(.385) 

.780 
(.238) 

18.053 
(.2  75) 

.800 
(.198) 

18.621 
(.255) 

.770 
(.265) 

t 

0.999 
L9.94 

5,536,000 
173 

20.152 
(.146) 

.950 
(.211) 

.920 
(.156) 

19.511 

(.138) 

.920 
(.252) 

19.979 
(.128) 

.950 

(.142) 
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